function [Pred_Marriage,TOTAL_SURPLUS,ALLSTORE,gamma,Surplus,Surplus_]=solveMatch20191030(Param,Glob,iprov,Period,Num,PolicyParam)


A=Param.A;
AgeParam=Param.AgeParam;
shadow=Glob.shadow;
shadow=shadow(iprov);
beta=Param.beta;
auxbeta=zeros(3,1);
auxbeta(1)=beta(1)+abs(beta(2))*shadow; % penalty for no access
auxbeta(2)=beta(3)+beta(4)*shadow; % payoff if access, men
auxbeta(3)=beta(5)+beta(6)*shadow; % payoff if access, women
lambdaAgeM=Param.lambdaAgeM;
lambdaEducF=Param.lambdaEducF;
lambdaH=Param.lambdaH;
lambdaHNI=Param.lambdaHNI;
sigma=Param.sigma;
delta=Param.delta;
OutsideCoeff=Param.Outside;
alpha=Param.alpha;
Nationality_=Param.Nationality;
lambdaNation=Param.lambdaNation;
OutsideNI=Param.OutsideNI;
if nargin==5
    PolicyParam=Glob.PolicyParam;
end
Policy=Glob.Policy;

seed=Glob.seed;
rng(seed+iprov+100*Num,'twister');
dimS=Glob.dimS;

dimAge=Glob.dimAge;
dimEduc=Glob.dimEduc;
dimNation=Glob.dimNation;
dimHouse=Glob.dimHouse;
ageGrid=Glob.ageGrid;
educGrid=Glob.educGrid;
ReducSize1=Glob.ReducSize1;
ReducSize2=Glob.ReducSize2;
LowParam=Glob.LowParam;
VecTotSingles=Glob.VecTotSingles;

if Period==1
DistMen=Glob.Singles_Men_Before;
DistMen=DistMen(:,iprov);
DistWomen=Glob.Singles_Women_Before;
DistWomen=DistWomen(:,iprov);
RatioSingleMF=Glob.RatioSingleMF_before;
RatioSingleMF=RatioSingleMF(iprov);RatioSingleMF=1;
OutCoeffGen=OutsideCoeff(iprov,1);
ProbHouse=Glob.ProbHouse;
ProbHouse=ProbHouse(:,[1 2]);
ShareSingle=VecTotSingles(:,iprov,1);
ShareSingle(1)=0;
elseif Period==2
DistMen=Glob.Singles_Men_After;
DistMen=DistMen(:,iprov);
DistWomen=Glob.Singles_Women_After;
DistWomen=DistWomen(:,iprov);
RatioSingleMF=Glob.RatioSingleMF_after;
RatioSingleMF=RatioSingleMF(iprov);RatioSingleMF=1;
OutCoeffGen=OutsideCoeff(iprov,2);
ProbHouse=Glob.ProbHouse;
ProbHouse=ProbHouse(:,[3 4]);
ShareSingle=VecTotSingles(:,iprov,2);
ShareSingle(1)=0;
Nationality_(1,3)=Nationality_(1,3)+lambdaNation(1);
Nationality_(3,1)=Nationality_(3,1)+lambdaNation(1);
end

if Policy==2  % Keep before initial distributions, but let beta vary depending on the period
DistMen=Glob.Singles_Men_Before;
DistMen=DistMen(:,iprov);
DistWomen=Glob.Singles_Women_Before;
DistWomen=DistWomen(:,iprov);
RatioSingleMF=Glob.RatioSingleMF_before;
RatioSingleMF=RatioSingleMF(iprov);RatioSingleMF=1;
OutCoeffGen=OutsideCoeff(iprov,1);
ProbHouse=Glob.ProbHouse;
ProbHouse=ProbHouse(:,[1 2]);
ShareSingle=VecTotSingles(:,iprov,1);
ShareSingle(1)=0;
Nationality_=Param.Nationality;
end

% dimS is the number of men, dimSW of women
dimSW=min(floor(dimS/RatioSingleMF),dimS);

educGrid2=kron(educGrid,ones(dimAge*dimNation,1));
educGrid3=kron(educGrid2,ones(dimHouse,1));
ageGrid2=kron(kron(ones(dimEduc,1),ageGrid),ones(dimNation,1));
ageGrid3=kron(ageGrid2,ones(dimHouse,1));
nationGrid2=kron(ones(dimAge*dimEduc,1),(1:dimNation)');
nationGrid3=kron(nationGrid2,ones(dimHouse,1));
HouseGrid3=kron(ones(dimAge*dimEduc*dimNation,1),(1:dimHouse)');

N=sum(DistMen.*(nationGrid2==1));
F=sum(DistMen.*(nationGrid2~=1));
gamma=alpha*N/(F+N);

%% Drawing random shocks

love=exp(randn(dimS,dimS)*sigma);
love_=exp(randn(dimS,dimS)*sigma);
uM=rand(dimS,1);
uF=rand(dimSW,1);
uM_=rand(dimS,1);
uF_=rand(dimSW,1);

rand1=exp(randn(1,dimS));
%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Pool of Italians without contact with foreigners
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

dimS=Glob.dimS;

auxItalianMen=DistMen.*(nationGrid2==1)*alpha;
auxItalianMen=[auxItalianMen.*(1-ProbHouse(:,1)) auxItalianMen.*ProbHouse(:,1)];
auxItalianMen=reshape(auxItalianMen',dimEduc*dimAge*dimNation*dimHouse,1);
auxDistMen=cumsum(auxItalianMen)./sum(auxItalianMen);

auxItalianWomen=DistWomen.*(nationGrid2==1)*alpha;
auxItalianWomen=[auxItalianWomen.*(1-ProbHouse(:,2)) auxItalianWomen.*ProbHouse(:,2)];
auxItalianWomen=reshape(auxItalianWomen',dimEduc*dimAge*dimNation*dimHouse,1);
auxDistWomen=cumsum(auxItalianWomen)./sum(auxItalianWomen);



auxtypeM=sortrows(uM,1);
auxtypeF=sortrows(uF,1);

typeM=sum(repmat(auxtypeM,[1 dimEduc*dimAge*dimNation*dimHouse])>repmat(auxDistMen',[dimS 1]),2)+1;
typeF=sum(repmat(auxtypeF,[1 dimEduc*dimAge*dimNation*dimHouse])>repmat(auxDistWomen',[dimSW 1]),2)+1;
typeF=[typeF;ones(dimS-dimSW,1)];

if Policy==6&&nationGrid3(PolicyParam(1))==1   % Adding an italian male to the first market
    love=[[love;rand1] exp(-10*ones(dimS+1,1))];
    dimS=dimS+1;
    dimSW=dimS;
    typeM(dimS,1)=PolicyParam(1);  % additional man
    typeF(dimS,1)=1; % does not matter. 
end

if Policy==6&&nationGrid3(PolicyParam(1))~=1   % Irrelevant guy just to make it square
    love=[[love;exp(-10*ones(1,dimS))] exp(-10*ones(dimS+1,1))];
    love(end,end)=0; % the irrelevant man and woman will match for sure see below
    dimS=dimS+1;
    dimSW=dimS;
    typeM(dimS,1)=1;  % does not matter
    typeF(dimS,1)=1; % does not matter
end


%% Age and education
X=[ageGrid3(typeF) educGrid3(typeF) HouseGrid3(typeF)];
X(dimSW+1:end,:)=0;
Y=[ageGrid3(typeM) educGrid3(typeM) HouseGrid3(typeM)];


indAgeM=sum(repmat(Y(:,1),[1 dimAge])>=repmat(ageGrid',[dimS 1]),2);
indAgeF=sum(repmat(X(:,1),[1 dimAge])>=repmat(ageGrid',[dimS 1]),2);
ind=indAgeM+(indAgeF'-1)*dimAge;
AgeEffect=AgeParam(ind);
%diffEduc=(Y(:,2)~=X(:,2)');
LowLow=(Y(:,2)==10)&(X(:,2)'==10);
LowHigh=(Y(:,2)==10)&(X(:,2)'==15);
HighLow=(Y(:,2)==15)&(X(:,2)'==10);
HighHigh=(Y(:,2)==15)&(X(:,2)'==15);
auxEduc=exp(LowLow*A(1)+LowHigh*A(2)+HighLow*A(3)+HighHigh*A(4));
LowHighH=(Y(:,3)==1)&(X(:,3)'==2);
HighLowH=(Y(:,3)==2)&(X(:,3)'==1);
HighHighH=(Y(:,3)==2)&(X(:,3)'==2);
HouseEffect=exp(LowHighH*lambdaH(1)+HighLowH*lambdaH(2)+HighHighH*lambdaH(3));


S=auxEduc.*AgeEffect.*HouseEffect;


%% Outside option

OutsideOptionM=repmat(OutCoeffGen,[1 dimS]);

OutsideOptionF=repmat(OutCoeffGen,[dimS 1]);

OutsideOption=OutsideOptionM/2+OutsideOptionF/2;


%% Global surplus
Surplus=S.*love-OutsideOption;
auxSurplus=S.*love;

NetSurplus=max(Surplus,0);
if Policy==6&&nationGrid3(PolicyParam(1))==1
    Surplus(:,dimS)=-100; % She remains single
    NetSurplus(:,dimS)=-100; % She remains single
end
if Policy==6&&nationGrid3(PolicyParam(1))~=1
    Surplus(:,dimS)=-100; % no one matches with the irrelevant woman
    Surplus(dimS,:)=-100; % no one matches with the irrelevant man
    Surplus(dimS,dimS)=-1; % the irrelevant man and woman match for sure but remain single
    NetSurplus(:,dimS)=-100; % no one matches with the irrelevant woman
    NetSurplus(dimS,:)=-100; % no one matches with the irrelevant man
    NetSurplus(dimS,dimS)=-1; % the irrelevant man and woman match for sure but remain single
end



matchresult=gurobSolve(NetSurplus,dimS);
match=reshape(matchresult.x,dimS,dimS);
[i1,OptWomen]=max(match,[],2);

ind=sub2ind([dimS dimS],(1:dimS)',OptWomen);
OptSurplus=Surplus(ind);
Singles=OptSurplus<0;
OptLove=auxSurplus(ind);

M=[typeM typeF(OptWomen) ageGrid3(typeM) ageGrid3(typeF(OptWomen)) educGrid3(typeM) educGrid3(typeF(OptWomen)) Singles OptSurplus];
totalSurplus=sum(auxSurplus.*match,2).*(1-Singles)+Singles.*sum(OutsideOption.*match,2);


%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% %% Italians and Foreigners
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

dimS=Glob.dimS;
dimS=floor(dimS*ReducSize2);
dimSW=min(floor(dimS/RatioSingleMF),dimS);

auxMen=DistMen.*(nationGrid2==1)*(1-alpha)+DistMen.*(nationGrid2~=1);
auxMen=[auxMen.*(1-ProbHouse(:,1)) auxMen.*ProbHouse(:,1)];
auxMen=reshape(auxMen',dimEduc*dimAge*dimNation*dimHouse,1);
auxDistMen=cumsum(auxMen)./sum(auxMen);

auxWomen=DistWomen.*(nationGrid2==1)*(1-alpha)+DistWomen.*(nationGrid2~=1);
auxWomen=[auxWomen.*(1-ProbHouse(:,1)) auxWomen.*ProbHouse(:,1)];
auxWomen=reshape(auxWomen',dimEduc*dimAge*dimNation*dimHouse,1);
auxDistWomen=cumsum(auxWomen)./sum(auxWomen);




auxtypeM_=sortrows(uM_,1);
auxtypeF_=sortrows(uF_,1);

typeM_=sum(repmat(auxtypeM_,[1 dimEduc*dimAge*dimNation*dimHouse])>repmat(auxDistMen',[dimS 1]),2)+1;
typeF_=sum(repmat(auxtypeF_,[1 dimEduc*dimAge*dimNation*dimHouse])>repmat(auxDistWomen',[dimSW 1]),2)+1;
typeF_=[typeF_;ones(dimS-dimSW,1)];


if Policy==1
    typeM_=[typeM_;ones(dimS*PolicyParam,1)*13];
    typeF_=[typeF_;ones(dimS*PolicyParam,1)*13];
    dimS=dimS*(1+PolicyParam);
    dimSW=dimSW*(1+PolicyParam);
end



if Policy==6   % Adding a male to the second market
    love_=[[love_;rand1] exp(-10*ones(dimS+1,1))];
    dimS=dimS+1;
    dimSW=dimS;
    typeM_(dimS,1)=PolicyParam(1);  
    typeF_(dimS,1)=1; 

end



if Policy==4 % Influx of young Africans 
    dimAfrican=PolicyParam(1);
    typeM_=[typeM_;ones(dimAfrican,1).*PolicyParam(2)];
    typeF_=[typeF_;ones(dimAfrican,1).*PolicyParam(2)];
    dimS2=dimS+dimAfrican;
    dimSW2=dimSW+dimAfrican;
    auxlove_=exp(randn(dimS2,dimS2)*sigma);
    auxlove_(1:dimS,1:dimS)=love_;
    love_=auxlove_;
    dimS=dimS2;
    dimSW=dimSW2;
end

if Policy==5 % Influx of Migrants
    dimNewMigrant=PolicyParam(1);
    iforeign=nationGrid2~=1&ageGrid2<=27&educGrid2==10; % Non Italians of age between 22 and 32, low educ
    auxDistMen=mean(DistMen(iforeign,:),2);
    auxNation=nationGrid2(iforeign);
    auxDistWomen=mean(DistWomen(iforeign,:),2);

    for i=1:7
    ind=auxNation==i+1;
    shareMigrantM(i)=sum(auxDistMen(ind))/sum(auxDistMen);
    shareMigrantF(i)=sum(auxDistWomen(ind))/sum(auxDistWomen);
    end
    CumShareMigrantM=cumsum(shareMigrantM);
    uOriginM=rand(dimNewMigrant,1);
    OriginMigrant=sum(repmat(uOriginM,[1 7])>repmat(CumShareMigrantM,[dimNewMigrant 1]),2)+2;
    OriginMigrant=OriginMigrant*0+PolicyParam(3);
    TypeSurgeM1=repmat(OriginMigrant,[1 dimAge*dimNation*dimEduc*dimHouse])==repmat(nationGrid3',[dimNewMigrant 1]);
    TypeSurgeM2=10*ones(dimNewMigrant,dimAge*dimNation*dimEduc*dimHouse)==repmat(educGrid3',[dimNewMigrant 1]);
    TypeSurgeM3=PolicyParam(2)*ones(dimNewMigrant,dimAge*dimNation*dimEduc*dimHouse)==repmat(ageGrid3',[dimNewMigrant 1]);
    TypeSurgeM4=ones(dimNewMigrant,dimAge*dimNation*dimEduc*dimHouse)==repmat(HouseGrid3',[dimNewMigrant 1]);
    TypeSurgeM=(TypeSurgeM1==1)&(TypeSurgeM2==1)&(TypeSurgeM3==1)&(TypeSurgeM4==1);
    [aux,TypeSurgeM]=max(TypeSurgeM,[],2);
    
    CumShareMigrantF=cumsum(shareMigrantF);
    uOriginF=rand(dimNewMigrant,1);
    OriginMigrant=sum(repmat(uOriginF,[1 7])>repmat(CumShareMigrantF,[dimNewMigrant 1]),2)+2;
    OriginMigrant=OriginMigrant*0+PolicyParam(3);
    TypeSurgeF1=repmat(OriginMigrant,[1 dimAge*dimNation*dimEduc*dimHouse])==repmat(nationGrid3',[dimNewMigrant 1]);
    TypeSurgeF2=10*ones(dimNewMigrant,dimAge*dimNation*dimEduc*dimHouse)==repmat(educGrid3',[dimNewMigrant 1]);
    TypeSurgeF3=PolicyParam(2)*ones(dimNewMigrant,dimAge*dimNation*dimEduc*dimHouse)==repmat(ageGrid3',[dimNewMigrant 1]);
    TypeSurgeF4=ones(dimNewMigrant,dimAge*dimNation*dimEduc*dimHouse)==repmat(HouseGrid3',[dimNewMigrant 1]);
    TypeSurgeF=(TypeSurgeF1==1)&(TypeSurgeF2==1)&(TypeSurgeF3==1)&(TypeSurgeF4==1);
    [aux,TypeSurgeF]=max(TypeSurgeF,[],2);
   
    typeM_=[typeM_;TypeSurgeM];
    typeF_=[typeF_;TypeSurgeF];
    dimS2=dimS+dimNewMigrant;
    dimSW2=dimSW+dimNewMigrant;
    auxlove_=exp(randn(dimS2,dimS2)*sigma);
    auxlove_(1:dimS,1:dimS)=love_;
    love_=auxlove_;
    dimS=dimS2;
    dimSW=dimSW2;
end



%% Legal access
NoAccessM_=nationGrid3(typeM_)>Period+1;
NoAccessF_=nationGrid3(typeF_)>Period+1;
aux=logical(NoAccessM_.*NoAccessF_');
Beta=ones(dimS,dimS);
Beta(aux)=auxbeta(1);
auxM=nationGrid3(typeM_)<=Period+1;
auxF=nationGrid3(typeF_)>Period+1;
aux=logical(auxM.*auxF');
Beta(aux)=auxbeta(2);
auxM=nationGrid3(typeM_)>Period+1;
auxF=nationGrid3(typeF_)<=Period+1;
aux=logical(auxM.*auxF');
Beta(aux)=auxbeta(3);


%% Age and education
X_=[ageGrid3(typeF_) educGrid3(typeF_) nationGrid3(typeF_) HouseGrid3(typeF_)];  
X_(dimSW+1:end,:)=0;

Y_=[ageGrid3(typeM_) educGrid3(typeM_)  nationGrid3(typeM_) HouseGrid3(typeM_)];


indAgeM=sum(repmat(Y_(:,1),[1 6])>=repmat(ageGrid',[dimS 1]),2);
indAgeF=sum(repmat(X_(:,1),[1 6])>=repmat(ageGrid',[dimS 1]),2);
ind=indAgeM+(indAgeF'-1)*dimAge;
AgeEffect_=AgeParam(ind);


LowLow=(Y_(:,2)==10)&(X_(:,2)'==10);
LowHigh=(Y_(:,2)==10)&(X_(:,2)'==15);
HighLow=(Y_(:,2)==15)&(X_(:,2)'==10);
HighHigh=(Y_(:,2)==15)&(X_(:,2)'==15);
auxEduc_=exp(LowLow*A(1)+LowHigh*A(2)+HighLow*A(3)+HighHigh*A(4));


%% Culture and language
nationM_=nationGrid3(typeM_);
nationF_=nationGrid3(typeF_);

indCult_=nationM_+(nationF_'-1)*dimNation;
NationalityDist_=Nationality_(indCult_);
NationalityAge_Educ=(ageGrid3(typeM_)-22).*(nationGrid3(typeM_)<=2).*(nationGrid3(typeF_)'>2)*lambdaAgeM... 
        +lambdaEducF*(ageGrid3(typeM_)-22).*(educGrid3(typeM_)==10).*(nationGrid3(typeM_)<=2).*(nationGrid3(typeF_)'>2);
        
NationalityDist_=NationalityDist_+NationalityAge_Educ;

lambda_=reshape(NationalityDist_,dimS,dimS);
lambda_=max(lambda_,0);

%% Housing
bothItalian=(Y_(:,3)<=2)&(X_(:,3)'<=2);
LowHighHS=(Y_(:,4)==1)&(X_(:,4)'==2)&not(bothItalian);
HighLowHS=(Y_(:,4)==2)&(X_(:,4)'==1)&not(bothItalian);
HighHighHS=(Y_(:,4)==2)&(X_(:,4)'==2)&not(bothItalian);
LowHighHII=(Y_(:,4)==1)&(X_(:,4)'==2)&bothItalian;
HighLowHII=(Y_(:,4)==2)&(X_(:,4)'==1)&bothItalian;
HighHighHII=(Y_(:,4)==2)&(X_(:,4)'==2)&bothItalian;

HouseEffect_=exp(LowHighHS*lambdaHNI(1)+HighLowHS*lambdaHNI(2)+HighHighHS*lambdaHNI(3)+LowHighHII*lambdaH(1)+HighLowHII*lambdaH(2)+HighHighHII*lambdaH(3));

S_=AgeEffect_.*auxEduc_.*HouseEffect_;


%% Outside option
OutsideOption_=OutCoeffGen+OutsideNI(nationM_)+OutsideNI(nationF_)' ...
            +lambdaNation(2)*ShareSingle(nationM_)+lambdaNation(2)*ShareSingle(nationF_)';
OutsideOptionM_=OutCoeffGen/2+OutsideNI(nationM_)+lambdaNation(2)*ShareSingle(nationM_);
OutsideOptionF_=OutCoeffGen/2+OutsideNI(nationF_)+lambdaNation(2)*ShareSingle(nationF_);
%% Global surplus
Surplus_=S_.*Beta.*lambda_.*love_-OutsideOption_;
auxSurplus_=S_.*Beta.*lambda_.*love_;

% Net surplus some small fixed noise so that the allocation is the same
% when the surplus is negative
NetSurplus_=max(Surplus_,0);
if Policy==6
    Surplus_(:,dimS)=-100; % She remains single
    NetSurplus_(:,dimS)=-100; % She remains single
end


matchresult_=gurobSolve(NetSurplus_,dimS);
match_=reshape(matchresult_.x,dimS,dimS);
[i1,OptWomen_]=max(match_,[],2);

ind=sub2ind([dimS dimS],(1:dimS)',OptWomen_);
OptSurplus_=Surplus_(ind);
SinglesM_=OptSurplus_<0;
% OptLove_=love_(ind);
OptLove_=auxSurplus_(ind);

M_=[typeM_ typeF_(OptWomen_) ageGrid3(typeM_) ageGrid3(typeF_(OptWomen_)) educGrid3(typeM_) educGrid3(typeF_(OptWomen_)) nationGrid3(typeM_) nationGrid3(typeF_(OptWomen_))...
    HouseGrid3(typeM_) HouseGrid3(typeF_(OptWomen_)) SinglesM_ OptSurplus_ OptWomen_];
totalSurplus_=sum(auxSurplus_.*match_,2).*(1-SinglesM_)+SinglesM_.*(OutsideOptionM_+OutsideOptionF_); 

Pred_Marriage__=zeros(dimEduc*dimAge*dimNation*dimHouse,dimEduc*dimAge*dimNation*dimHouse);
Pred_Marriage_=zeros(dimEduc*dimAge*dimNation*dimHouse,dimEduc*dimAge*dimNation*dimHouse);
m=kron(eye(dimEduc*dimAge*dimNation),ones(dimHouse,1));
% Marriages Italians only 
Singles=OptSurplus<LowParam;
q=sub2ind([dimEduc*dimAge*dimNation*dimHouse dimEduc*dimAge*dimNation*dimHouse],M(not(Singles),1),M(not(Singles),2));
tab=accumarray(q,1);
Pred_Marriage__(1:length(tab))=tab;
Pred_Marriage__=m'*Pred_Marriage__*m/(dimS*ReducSize1);

% Marriages Italians + foreigners
SinglesM_=OptSurplus_<LowParam;
q=sub2ind([dimEduc*dimAge*dimNation*dimHouse dimEduc*dimAge*dimNation*dimHouse],M_(not(SinglesM_),1),M_(not(SinglesM_),2));
tab_=accumarray(q,1);
Pred_Marriage_(1:length(tab_))=tab_;
Pred_Marriage_=m'*Pred_Marriage_*m/(dimS*ReducSize2);


Pred_Marriage=Pred_Marriage__*gamma+(1-gamma)*Pred_Marriage_;
TOTAL_SURPLUS=[totalSurplus;totalSurplus_];
STORE=[nationGrid3(typeM) nationGrid3(typeF(OptWomen)) ageGrid3(typeM) ageGrid3(typeF(OptWomen)) educGrid3(typeM) educGrid3(typeF(OptWomen))... 
Singles OptWomen OutsideOptionM'/2 OutsideOptionF/2 sum(auxSurplus.*match,2).*(1-Singles) typeM typeF(OptWomen) OptLove];

STORE_=[nationGrid3(typeM_) nationGrid3(typeF_(OptWomen_)) ageGrid3(typeM_) ageGrid3(typeF_(OptWomen_)) educGrid3(typeM_) educGrid3(typeF_(OptWomen_))...
    SinglesM_ OptWomen_ OutsideOptionM_ OutsideOptionF_(OptWomen_) sum(auxSurplus_.*match_,2).*(1-SinglesM_) typeM_ typeF_(OptWomen_) OptLove_];


ALLSTORE.STORE=STORE;
ALLSTORE.STORE_=STORE_;
STORE3YC=NaN;
STORE3YC_=NaN;
STORE3YB=NaN;
STORE3YB_=NaN;
STORE3YA=NaN;
STORE3YA_=NaN;

%% 3 Years after ...
if Period==1&Policy==0
    % natives before and after
    REL=rand(dimS,1)>delta(3);
    ind=sub2ind([dimS dimS],(1:dimS)',OptWomen);
    indd=ind(Surplus(ind)>0);
    love3=love;
    innovation=randn(dimS,dimS)*delta(1);
    love3(indd)=love3(indd)+innovation(indd);
    Surplus3=S.*love3-OutsideOption;
    Divorce3YB=(Surplus3(ind)<0)&(Singles==0)&(REL==0);
    Married3YB=Singles==0;
    STORE3YB=[nationGrid3(typeM) nationGrid3(typeF(OptWomen)) ageGrid3(typeM) ageGrid3(typeF(OptWomen)) educGrid3(typeM) educGrid3(typeF(OptWomen)) Divorce3YB Married3YB];
    STORE3YA=STORE3YB;
    
    %  foreigners before
    REL_=rand(dimS,1)>delta(4);
    ind=sub2ind([dimS dimS],(1:dimS)',OptWomen_);
    indd=ind(Surplus_(ind)>0);
    RAND=randn(dimS,dimS);
    love3_=love_;
    innovation_=RAND*delta(2);
    love3_(indd)=love3_(indd)+innovation_(indd);
    Surplus3_=S_.*Beta.*lambda_.*love3_-OutsideOption_;
    Divorce3YB_=(Surplus3_(ind)<0)&(SinglesM_==0)&(REL_==0);    
    Married3YB_=SinglesM_==0;
    
    STORE3YB_=[nationGrid3(typeM_) nationGrid3(typeF(OptWomen_)) ageGrid3(typeM_) ageGrid3(typeF(OptWomen_)) educGrid3(typeM_) educGrid3(typeF(OptWomen_)) Divorce3YB_ Married3YB_];
    
    %  foreigners after
    NoAccessM_=nationGrid3(typeM_)>3; % This cohort experiences the enlargement 
    NoAccessF_=nationGrid3(typeF_)>3;
    aux=logical(NoAccessM_.*NoAccessF_');
    BetaA=ones(dimS,dimS);
    BetaA(aux)=auxbeta(1);
    auxM=nationGrid3(typeM_)<=3;
    auxF=nationGrid3(typeF_)>3;
    aux=logical(auxM.*auxF');
    BetaA(aux)=auxbeta(2);
    auxM=nationGrid3(typeM_)>3;
    auxF=nationGrid3(typeF_)<=3;
    aux=logical(auxM.*auxF');
    BetaA(aux)=auxbeta(3);
    
    Surplus3A_=S_.*BetaA.*lambda_.*love3_-OutsideOption_;
    Divorce3YA_=(Surplus3A_(ind)<0)&(SinglesM_==0)&(REL_==0);    
    Married3YA_=SinglesM_==0;
    STORE3YA_=[nationGrid3(typeM_) nationGrid3(typeF(OptWomen_)) ageGrid3(typeM_) ageGrid3(typeF(OptWomen_)) educGrid3(typeM_) educGrid3(typeF(OptWomen_)) Divorce3YA_ Married3YA_];
end
if Period==2
    % natives
    REL=rand(dimS,1)>delta(3);
    ind=sub2ind([dimS dimS],(1:dimS)',OptWomen);
    indd=ind(Surplus(ind)>0);
    love3=love;
    innovation=randn(dimS,dimS)*delta(1);
    love3(indd)=love3(indd)+innovation(indd);
    Surplus3=S.*love3-OutsideOption;
    Divorce3YC=(Surplus3(ind)<0)&(Singles==0)&(REL==0);
    Married3YC=Singles==0;
    STORE3YC=[nationGrid3(typeM) nationGrid3(typeF(OptWomen)) ageGrid3(typeM) ageGrid3(typeF(OptWomen)) educGrid3(typeM) educGrid3(typeF(OptWomen)) Divorce3YC Married3YC];
       
    
    %  foreigners 
    NoAccessM_=nationGrid3(typeM_)>3; % This cohort experiences the enlargement 
    NoAccessF_=nationGrid3(typeF_)>3;
    aux=logical(NoAccessM_.*NoAccessF_');
    BetaA=ones(dimS,dimS);
    BetaA(aux)=auxbeta(1);
    auxM=nationGrid3(typeM_)<=3;
    auxF=nationGrid3(typeF_)>3;
    aux=logical(auxM.*auxF');
    BetaA(aux)=auxbeta(2);
    auxM=nationGrid3(typeM_)>3;
    auxF=nationGrid3(typeF_)<=3;
    aux=logical(auxM.*auxF');
    BetaA(aux)=auxbeta(3);
    
    REL_=rand(dimS,1)>delta(4);
    ind=sub2ind([dimS dimS],(1:dimS)',OptWomen_);
    indd=ind(Surplus_(ind)>0);
    RAND=randn(dimS,dimS);
    love3_=love_;
    innovation_=RAND*delta(2);
    love3_(indd)=love3_(indd)+innovation_(indd);

    Surplus3C_=S_.*BetaA.*lambda_.*love3_-OutsideOption_;
    Divorce3YC_=(Surplus3C_(ind)<0)&(SinglesM_==0)&(REL_==0);    
    Married3YC_=SinglesM_==0;
    STORE3YC_=[nationGrid3(typeM_) nationGrid3(typeF(OptWomen_)) ageGrid3(typeM_) ageGrid3(typeF(OptWomen_)) educGrid3(typeM_) educGrid3(typeF(OptWomen_)) Divorce3YC_ Married3YC_];

end
ALLSTORE.STORE3YC=STORE3YC;   
ALLSTORE.STORE3YC_=STORE3YC_;
ALLSTORE.STORE3YB=STORE3YB;
ALLSTORE.STORE3YB_=STORE3YB_;
ALLSTORE.STORE3YA=STORE3YA;
ALLSTORE.STORE3YA_=STORE3YA_;


end


